function getImpulseResponses_yieldCurveFactors(outEstim,IRFlength,shockSizeScalar,xfGIRF,xsGIRF)
%% Impulse response functions
model       = outEstim.model;
ne          = size(model.eta,2);
vectorMom3  = model.vectorMom3;
reduceStatesON = 1;
selectY = [1 2 3 5];
m = outEstim.setupFilter.CSregress_m;
[model,errorMes]  = perturbationDSGEmodel(model.params,model.orderApp,outEstim.setupFilter,selectY,reduceStatesON);
% Level factor: (r^(1) + r^(8) + r^(40))/3
model.labelYFactor        = [{'$Level_t$'}];
model.Yfactor_x(1,:)      = 4*(model.byx(1,:)  + model.byx(8,:)  + model.byx(40,:))/3;
model.Yfactor_xx(1,:,:)   = 4*reshape(model.byxx(1,:)  + model.byxx(8,:)  + model.byxx(40,:),model.nx,model.nx)/3;
model.Yfactor_xxx(1,:,:,:)= 4*reshape(model.byxxx(1,:) + model.byxxx(8,:) + model.byxxx(40,:),model.nx,model.nx,model.nx)/3;
model.Yfactor_ss(1,:)     = 4*(model.byss(1,:)  + model.byss(8,:)  + model.byss(40,:))/3;
model.Yfactor_ssx(1,:)    = 4*(model.byssx(1,:)  + model.byssx(8,:)  + model.byssx(40,:))/3;
model.Yfactor_sss(1,:)    = 4*(model.bysss(1,:)  + model.bysss(8,:)  + model.bysss(40,:))/3;

% Slope factor: r^(40) - r^(1) 
model.labelYFactor        = [model.labelYFactor, {'$Slope_t$'}];
model.Yfactor_x(2,:)      = 4*(model.byx(40,:)  - model.byx(1,:))/3;
model.Yfactor_xx(2,:,:)   = 4*reshape(model.byxx(40,:)  - model.byxx(1,:),model.nx,model.nx)/3;
model.Yfactor_xxx(2,:,:,:)= 4*reshape(model.byxxx(40,:) - model.byxxx(1,:),model.nx,model.nx,model.nx)/3;
model.Yfactor_ss(2,:)     = 4*(model.byss(40,:)  - model.byss(1,:))/3;
model.Yfactor_ssx(2,:)    = 4*(model.byssx(40,:) - model.byssx(1,:))/3;
model.Yfactor_sss(2,:)    = 4*(model.bysss(40,:) - model.bysss(1,:))/3;

% Curvature factor: r^(8)-r^(1) - (r^(40)-r^(8)) = 2*r^(8) - r^(1) - r^(40)
model.labelYFactor        = [model.labelYFactor, {'$Curvature_t$'}];
model.Yfactor_x(3,:)      = 4*(2*model.byx(8,:) - model.byx(1,:) - model.byx(40,:))/3;
model.Yfactor_xx(3,:,:)   = 4*reshape(2*model.byxx(8,:)  - model.byxx(1,:) - model.byxx(40,:),model.nx,model.nx)/3;
model.Yfactor_xxx(3,:,:,:)= 4*reshape(2*model.byxxx(8,:) - model.byxxx(1,:) - model.byxxx(40,:),model.nx,model.nx,model.nx)/3;
model.Yfactor_ss(3,:)     = 4*(2*model.byss(8,:)  - model.byss(1,:) - model.byss(40,:))/3;
model.Yfactor_ssx(3,:)    = 4*(2*model.byssx(8,:)  + model.byssx(1,:) - model.byssx(40,:))/3;
model.Yfactor_sss(3,:)    = 4*(2*model.bysss(8,:)  + model.bysss(1,:) - model.bysss(40,:))/3;

% Benchmark with GAMApai = 0
params0 = model.params;
params0.OMEGAd = 0;
[model0,errorMes]       = perturbationDSGEmodel(params0,model.orderApp,outEstim.setupFilter,selectY,reduceStatesON);
% Level factor: (r^(1) + r^(8) + r^(40))/3
model0.labelYFactor        = [{'$Level_t$'}];
model0.Yfactor_x(1,:)      = 4*(model0.byx(1,:)  + model0.byx(8,:)  + model0.byx(40,:))/3;
model0.Yfactor_xx(1,:,:)   = 4*reshape(model0.byxx(1,:)  + model0.byxx(8,:)  + model0.byxx(40,:),model0.nx,model0.nx)/3;
model0.Yfactor_xxx(1,:,:,:)= 4*reshape(model0.byxxx(1,:) + model0.byxxx(8,:) + model0.byxxx(40,:),model0.nx,model0.nx,model0.nx)/3;
model0.Yfactor_ss(1,:)     = 4*(model0.byss(1,:)  + model0.byss(8,:)  + model0.byss(40,:))/3;
model0.Yfactor_ssx(1,:)    = 4*(model0.byssx(1,:)  + model0.byssx(8,:)  + model0.byssx(40,:))/3;
model0.Yfactor_sss(1,:)    = 4*(model0.bysss(1,:)  + model0.bysss(8,:)  + model0.bysss(40,:))/3;

% Slope factor: r^(40) - r^(1) 
model0.labelYFactor        = [model0.labelYFactor, {'$Slope_t$'}];
model0.Yfactor_x(2,:)      = 4*(model0.byx(40,:)  - model0.byx(1,:))/3;
model0.Yfactor_xx(2,:,:)   = 4*reshape(model0.byxx(40,:)  - model0.byxx(1,:),model0.nx,model0.nx)/3;
model0.Yfactor_xxx(2,:,:,:)= 4*reshape(model0.byxxx(40,:) - model0.byxxx(1,:),model0.nx,model0.nx,model0.nx)/3;
model0.Yfactor_ss(2,:)     = 4*(model0.byss(40,:)  - model0.byss(1,:))/3;
model0.Yfactor_ssx(2,:)    = 4*(model0.byssx(40,:) - model0.byssx(1,:))/3;
model0.Yfactor_sss(2,:)    = 4*(model0.bysss(40,:) - model0.bysss(1,:))/3;

% Curvature factor: r^(8)-r^(1) - (r^(40)-r^(8)) = 2*r^(8) - r^(1) - r^(40)
model0.labelYFactor        = [model0.labelYFactor, {'$Curvature_t$'}];
model0.Yfactor_x(3,:)      = 4*(2*model0.byx(8,:) - model0.byx(1,:) - model0.byx(40,:))/3;
model0.Yfactor_xx(3,:,:)   = 4*reshape(2*model0.byxx(8,:)  - model0.byxx(1,:) - model0.byxx(40,:),model0.nx,model0.nx)/3;
model0.Yfactor_xxx(3,:,:,:)= 4*reshape(2*model0.byxxx(8,:) - model0.byxxx(1,:) - model0.byxxx(40,:),model0.nx,model0.nx,model0.nx)/3;
model0.Yfactor_ss(3,:)     = 4*(2*model0.byss(8,:)  - model0.byss(1,:) - model0.byss(40,:))/3;
model0.Yfactor_ssx(3,:)    = 4*(2*model0.byssx(8,:)  + model0.byssx(1,:) - model0.byssx(40,:))/3;
model0.Yfactor_sss(3,:)    = 4*(2*model0.bysss(8,:)  + model0.bysss(1,:) - model0.bysss(40,:))/3;

% Test the size of xfGIRF and xsGIRF, the point around which the GIRFs are
% computed
nx  = size(model.hx,1);
if ~isempty(xfGIRF)
    if size(xfGIRF,1) ~= nx
        error('xfGIRF must be of size nx*1')
    end
end
if ~isempty(xsGIRF)
    if size(xsGIRF,1) ~= nx
        error('xsGIRF must be of size nx*1')
    end
end

IRFyAll = zeros(3,IRFlength,3,ne);
for j=1:ne
    shockSize      = zeros(ne,1);
    shockSize(j,1) = shockSizeScalar;
    [IRFy,IRFx] = GIRFPruning_v2(model.Yfactor_x,model.Yfactor_xx,model.Yfactor_ss,model.Yfactor_xxx,model.Yfactor_ssx,model.Yfactor_sss,...
        model.hx,model.hxx,model.hss,model.hxxx,model.hssx,model.hsss,model.eta,...
        shockSize,vectorMom3,xfGIRF,xsGIRF,model.orderApp,IRFlength);
    
    % Benchmark: with GAMApai = 0
    [IRFy0,IRFx0] = GIRFPruning_v2(model0.Yfactor_x,model0.Yfactor_xx,model0.Yfactor_ss,model0.Yfactor_xxx,model0.Yfactor_ssx,model0.Yfactor_sss,...
        model0.hx,model0.hxx,model0.hss,model0.hxxx,model0.hssx,model0.hsss,model0.eta,...
        shockSize,vectorMom3,xfGIRF,xsGIRF,model.orderApp,IRFlength);
    
    % We store the Benchmark at pos 1 in IRFy and IRFx
    IRFy(:,:,1) = IRFy0(:,:,model.orderApp);
    
    IRFyAll(:,:,:,j) = IRFy;
end
model.labelne(1,2) = {'$d_t$'};
plotImpulseResponse_yieldCurveFactors(IRFyAll*100,model.labelYFactor,model.labelne,shockSize,model.orderApp);
end
